Semi-classical statistical approach to Frohlich condensation theory 

Jordane Preto 1 '!!] 

1 Aix-Marseille University, Campus de Luminy, case 907, CNRS Centre 
de Physique Theorique, UMR 7332, 13288 Marseille Cedex 09, France 

(Dated: February 25, 2013) 

Abstract 

O 

^vq Frohlich model equations describing phonon condensation in open systems of biological relevance 

are here reinvestigated in a semi-classical non-equilibrium statistical context (with "semi-classical" 
it is meant that the evolution of the system is described by means of classical equations with 
the addition of energy quantization). In particular, the assumptions that are necessary to deduce 
Frohlich rate equations are highlighted and we show how these hypotheses led us to write an ap- 
propriate form for the master equation. As a comparison with known previous results, analytical 
relations with the Wu- Austin quantum Hamiltonian description are emphasized. Finally, we show 
how solutions of the master equation can be implemented numerically and outline some representa- 
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tive results of the condensation effect. Our approach thus provides more information with respect 
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to the existing ones, in what we are concerned with the time evolution of the probability density 
functions instead of following average quantities. 
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I. INTRODUCTION 

Some decades ago, the study of open systems far from thermodynamic equilibrium 
showed, under suitable conditions, the emergence of self-organization. Striking similarities 
were observed among very different physical systems having in common the fact of being 
composed of many non-linearly interacting subsystems. When a control parameter, typi- 
cally the energy input rate, exceeds a critical value then the subsystems act cooperatively 
to self-organize in what is commonly referred to as a non-equilibrium phase transition. This 
is at variance with equilibrium phase transitions, for which the transition from disorder to 
order is driven by a lowering of temperature - the control parameter -below a critical value. 
Paradigmatic examples of non-equilibrium phase transitions are provided by the laser tran- 
sition and by the Rayleigh-Benard convective instability, a quantum and a classical system, 
respectively. In the early '80s of the last century, the emergence of collective properties 
in open systems was studied in a general and interdisciplinary framework referred to as 
"Synergetics" [lj. This fascinating topic was pioneered in the late '60s by H. Frohlich [2J. 

Frohlich considered a system consisting in z normal modes of frequency Uk (oj\ < 002 < 
... < u z ), each characterized by a discrete energy spectrum (phonons). Due to interactions 
with a surrounding heat bath, it was supposed that the dynamics of phonons in each mode 
was governed by several kinds of stochastic events : (i) linear events resulting in the absorp- 
tion or the emission of a single phonon from a particular normal mode, (ii) non-linear events 
resulting in the simultaneous absorption of one phonon from a normal mode, and emission 
of one phonon to another mode. Besides interactions with the heat bath, each normal mode 
is assumed to be fed energy by some external source. Denoting by (Nk) the average number 
of phonons in the mode of frequency Uk, Frohlich suggested that the dynamics of the system 
was described by the following rate equations 

^ = s k + MP) (W + 1 - (N k )e^) + 

(1) 
J> fci G0) [(N k + l){N d ) - (N k )(N j + l) e /^-^)] , k = 1, ..., z , 

where /3 = 1/kT with T the temperature of the heat bath. Here, the first r.h.s. Sk denotes 
the rate of external energy supply to the mode k whereas the second and the third terms 
account for the rate of change due to linear (i) and non-linear (ii) events, respectively (with 



appropriate temperature-dependent coupling constants <pk and Akj, respectively). Frohlich 
solved analytically Eq. (??) in the stationary case — - — = 0, showing that the total number 
of phonons increased linearly by increasing the total rate of external energy supply J^ fc S& 
beyond a certain threshold, while non-linear events (ii) tend to redistribute the energy excess 
into the mode of lowest frequency In other words, for sufficiently large value of the s^'s, the 
stationary state is shown to satisfy 

(Aq) ~ (N) 

z 

where (N) = ^(A/j), i.e., it can be considered that all the components of the system 
oscillate in a collective way at the lowest frequency of the spectrum. This state of average 
excitation is better known as Frohlich condensation. 

Since Frohlich's first proposal [2], Frohlich condensation has been investigated by many 
authors [3 J . In particular, a generic model of Hamiltonian was provided by Wu and Austin 
to derive Frohlich rate equations [Eqs. (??)] from a microscopic quantum basis [U Ej. This 
type of Hamiltonian was then found to induce excitations propagating coherently in a similar 
way as Davydov's solitons [6j [7] . From a practical point of view, Frohlich condensation has 
been suggested to play a central role in the self-organization of biological polar structures 
(here the set of normal modes arises from polar oscillations). This is the case, for instance, of 
microtubules and cell membranes that fulfilled the main requirements of Frohlich systems [3J 
[S] (there the hydrolysis of adenosine triphosphate (ATP) or guanosine triphosphate (GTP) 
represents a significant source of external energy supply). Recently, Pokorny experimentally 
observed a strong excitation in the spectrum of vibration of microtubules localized in the 
10 MHz range [9]. Among other applications of condensation, Frohlich suggested that the 
excitation of a particular mode of polar oscillations in macromolecular systems could lead to 
strong long-range dipole-dipole interactions. Applied to biological systems, it was expected 
that such forces would have a profound influence on the displacement of specific biological 
entities [TO] 111], and thus on the initiation of a particular cascade of chemical events. While 
long-range interactions have been reported at the cellular level [12], their possible role at 
the biomolecular level still remains an open question [13J. 

Now, coming back to Frohlich theory, it should be stressed that most studies performed 
on the condensation - including the derivation of Eqs. (??) - are based on a quantum 
Hamiltonian description akin to the one proposed originally by Wu and Austin. Given 



the semi-classical nature of the rate equations postulated by Frohlich, one can legitimately 
wonder how necessary a quantum description is to account for Frohlich condensation or even 
what are the ingredients needed to contemplate such a phenomenon. In the present paper, 
we propose to tackle this issue from a semi-classical point of view (here "semi-classical" 
means that the evolution of the system is described by means of classical equations with 
the addition of energy quantization). In section [III Frohlich rate equations (??) are derived 



from clear semi-classical arguments (II A and II B) and we show, in qualitative terms, how 



Frohlich condensation can be "deduced" from these arguments (II C). In section III, the 
classical master equation which is at the grounds of the rate equations, is given from the main 
results found in section |TTJ Then, the same equation is found to arise from the Wu- Austin 
Hamiltonian description under appropriate decoherence assumptions. Finally, as another 
original approach, Frohlich condensation is statistically studied by estimating the solution 
of the master equation numerically [section IV] . A lot of informations is thus provided that 
could, in our opinion, be of particular relevance to identify the phenomenon experimentally. 

II. DERIVATION OF THE RATE EQUATIONS 

For the sake of readability, the derivation of Frohlich rate equations has been itself splitted 
into two subsections. The first one (A) deals specifically with the derivation of the second 
r.h.s. of Eq. (??), i. e. that part which accounts for linear events, whereas the second one (B) 
is about the derivation of the equations in their entirety. In the former case, it is enough to 
consider the presence of a single normal mode in the system. Thus, this will provide a good 
introduction to the second subsection where the dynamics of all normal modes must be taken 
into account simultaneously in order to describe non-linear events. Finally, a subsection (C) 
has been added as a complement to show how Frohlich condensation can be "predicted" 
qualitatively on the basis of the transition probabilities given in (A) and (B) (see below). 

A. One variable process 

Although that point is usually not explicitly mentioned in the literature, Frohlich equa- 
tions are primarily based on the assumption that the time dependence of the number of 
phonons in each normal mode can be described as a homogeneous Markov process. Here 



one normal mode of frequency u k is considered and we call Nk(t) the stochastic process that 
accounts for the number of phonons in this mode. Under homogeneous Markov assumption 
it is possible to work with time-translation invariant conditional probabilities p(n k ,t\n' k ) 
such that 

p{n k ,t\n' k ) = Pioh{N k {t + t') = n k \N k (t') = n' k ,W > 0}, 
satisfies the following master equation [1] 



d t p(n k ,t\n k ) = ^ W{n k \m k )p{m k , t\n' k ) - W(m k \n k )p(n k ,t\n k ). (2) 

m k 

Here, the transition probabilities 

Tj rr/ | n ,• p{n k ,At\m k ) 
W{n k \m k ) = hm — , 

Ai->0 At 

are supposed to be well defined. Moreover, since Eq. (??) holds true irrespectively of the 
initial condition n' k , conditional probabilities p(...,t\n' k ) will be substituted with standard 
ones p(...,t) in what follows. 

As mentioned in the heading of the section, only interactions that lead to the absorption 
or the emission of one phonon at a time are considered here; thus W(n k \m k ) is zero except 
when m k = n k ±l or m k = n k {birth and death process). Then, the equation for the evolution 
of the average number of phonons in the mode k is easily deduced 

, k = ^ n k d t p(n k , t) 

n k >0 

= ^2 n k [W(n k \n k - l)p(n k - 1, t) + W(n k \n k + l)p{n k + 1, t) (3) 

ra fc >0 

-{W{n k + l\n k ) + W(n k - l\n k ))p(n kl t)\ . 
With some appropriate substitutions in the sum index with suitable relabeling, one gets 
^A = J2 [W(n k + IK) - W(n k - l\n k )]p(n k ,t), (4) 

n k >0 

where the boundary conditions 



W(-1|0) = 0, and p(n k ,t) = 0, if n k < 0, (5) 

have been used. 

Now, let us look for an appropriate expression for W(rik + l\n k ) and W(rik — l\n k ) in 
terms of n k . First, we suppose that all phonons are "independent" so that 

W(n k - l\n k ) = a k ((3)n k , (6) 

where ot k (f$) is the probability per time unit (which is temperature-dependent in general) 
that one phonon is emitted from the mode k. Second, W(n k + l\n k ) can be obtained by 
noticing that a birth and death process following boundary conditions (??) has a single 
stationary solution p s (n k ) through Eq. (??) (see the Appendix), and that this solution 
satisfies the detailed balance condition 

W{n k + l\n k )p s (n k ) = W(n k \n k + l)p s {n k + 1), V n k > 0. (7) 

Moreover, due to the contact between the system and the surrounding heat bath, it is 
required that the stationary solution is given by a Bose-Einstein distribution 

p s (n k ) = e - fihUknk /Z, where Z = ^ e~^ nk (8) 

ra fe >0 

is the one-normal mode partition function. Eq. (??) then leads to : 

W(n k + IK) = W(n k \n k + l) e "^* ( = ?) a k ((3)(n k + l)e~^K (9) 

Substituting this last equation and Eq. (??) into Eq. (??), we finally get the second 
r.h.s. of Eq. (??) : 

^ = MP) {(N k ) + 1 - (N k )e^) , (10) 

where we have let (p k (j3) = a k (p)e~ phWk . 

B. Many variables process 

We now want to consider a possible influence of all other normal modes on the mode of 
frequency u k . In this way, let N(t) = (Ni(t), N 2 (t), ..., N z (t)) be the vector of all Marko- 
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vian processes describing the number of phonons in each normal mode, so that the con- 
ditional probabilities are noted as p(n,t\n'), with n, n' G W, n = (ni, n 2 , ...,n z ) and 



n' = (n\,n' n' 



D "*2> •••■> lb zJ- 

In anticipation of a large number of interactions between the normal modes and the 
heat bath, we call R^, p, = l,...,s, all possible events that account for the simultaneous 
absorption(s) and/or emission(s) of phonon(s) from particular modes. Then, each event is 
characterized by a vector r M whose element rf (1 < i < z) is equal to the corresponding 
variation in the number of phonons in the mode of frequency Wj. Thus, the master equation 
takes the usual form 



dtp(n,t\n') = y^W^nln - r^)p(n - r M ,£|n') - W(n - r M |n)p(n,t|n / ). (11) 

In this way, the evolution of the average number of phonons in the mode k is immediately 
deduced (similarly to the previous subsection, initial conditions are omitted) 



d(N k ) 



^^2n k [W{n\n-r^)p(n-r^,t)-W{n-r^\n)p(n,t)}, (12) 

where the sum ^2 corresponds to J2 J2 ■■■ J2 ■ Besides, since events which do not 

n m>0ri2>0 n z >0 

result in a variation in the number of phonons in the mode k have zero contribution to the 
equation, Eq. (??) can be readily simplified. 

1. Second r.h.s. of Eq. (??) : as mentioned above, each event Rf 1 corresponds to the 
absorption or the emission of one phonon from a particular mode. Thus, there are only 
two events involving the mode k (2z events overall) and they are given by r^ = ±5jj~, 
for j = 1, ...,z. This case corresponds to the case tackled in the previous subsection 
whereby equation (??) was finally found. 

2. Third r.h.s. of Eq. (??) : here each event _R M corresponds to the simultaneous absorp- 
tion of one phonon from a normal mode, and emission of one phonon from another 
mode. In particular, with regard to the mode k, there are 2(z — 1) events (2z(z — 1) 
events overall), each characterized by r M such that 

k j 

I I (13) 

r^ = (0,....,0,±l,0,....,0, T l,0,...,0), Vj^h, je{i,zh 



where j G [1, z\ means j = 1, 2, ..., z. For the sake of clarity in what follows, we shall 
note W n (n k ;rij\ n k ± l;rij =F 1) the transition probabilities VF(n|n — r^) associated 
with these events. Likewise, p(n — r M , t) will be written as p n {n k ± 1; rij =F 1, t) (in this 
way, p(n,t) may be sometimes noted as p n (n k ;rij,t)). In this way, Eq. (??) becomes 

-j k = ^2^2 n k [W n (n k ;nj\n k - l;rij + l)p n (n k - I; rij + l,t) + 

W n (n k ;n j \n k + l;rij - l)p n (n k + l;rij - l,t) - 

W„(njfc + 1; rij - l\n k ; n j )p n (n k \ rij, t) - 

W n (n k - 1; nj + l|n fc ; nj)p n (n k ; rij, t)] . 

(14) 

Again, some appropriate substitutions on the sum indexes with suitable relabeling 
lead to 

d(N k ) 

W n (n k - 1; rij + l|n fc ; n,-)] p n (n fe ; rij, t), 



- ^2^2[W n (n k + \\rij - l\n k ;nj)- 

M rtk » (15) 



where the boundary conditions 

W n (—l;rij + l|0;rij) = 0, W n (n k + l;-l\n k ;0) = 0, and 

(16) 

Pn(n fc ; n i? t) = 0, if n fc or rij < 0, 

have been used. 

As previously, we look for the expression of the transition probabilities. On each 
normal mode j ^ k, we suppose that the phonons have the same probability 7fcj(/3, n k ) 
of being emitted while a phonon is absorbed in the mode k (in general terms, this 
quantity depends on the temperature as well as on the current number n k of phonons 
in the mode k). We thus get 

W n (n k + l;rij- l\n k ; rij) = W n (rij -l;n k + l\rij\ n k ) = -f kj {(3, n k )rij, 

(17) 

V n k , rij > and Vj^fe, jG [1,-j]. 
8 



At this stage, it is worth mentioning that the transition probabilities considered thus 
far have been given in such a way that the graph of the master equation is connected. 
Therefore, according to general results on master equations [13] , it can be found that 
a stationary solution of Eq (??) exists and is unique. Moreover, in so far as no other 
external source is present, a stationary solution for which detailed balanced condition 
is fulfilled, always exists. In our case, this solution is thus unique. Besides, since the 
system of normal modes interacts with the heat bath only, the stationary solution has 
to be given by a Bose-Einstein distribution. Thus 

W n (n k + 1; rij - l\n k ; nj)p s n (n k ] rij) = W n (n k ; 7ij\n k + 1; rij - l)p s n {n k + 1; rij - 1), 

V n k , rij > and V j ^ k, j € [1, z\ , 



(18) 



with p s n (n k - n 3 ) = J J e~^ ini /Z, where Z = \[J2 e 



i i rii>0 

Using Eq. (??), Eq. (??) becomes 



n k + 1 rij 

which must depend on k and j only (i.e., not on n k nor on rij). We finally deduce the 
properties 

JkjW, n k ) = A kj (P)(n k + 1), (19) 

with 

A jk ((3) = A kj ((3)e^-^\ Vj^k, je [1, z\ . (20) 

Using these equations and Eq. (??) in Eq. (??), we obtain 
d(N h ) 



dt m 



X)A fci (/3) [((N k + l)Nj) - (N k (N 3 + l)) e ^K-i)] , (21) 



which is exactly the third r.h.s. of Frohlich equations (??) in so far as correlations 
between the normal modes can be neglected, i.e. (N k (Nj + 1)) ~ (N k )(Nj + 1), or 
more specifically p(n, t) ~ ]Tp("-j,t). 



3. First r.h.s. of Eq. (??) (source term) : the constant rate of external energy supply 
may be worked out by considering the external source as a heat bath with infinite 
temperature, with which the system interacts linearly. Thus, the rate of change in the 
mode k due to the source is given by an equation similar to (??) 



^ = tf (A) (W + 1 - (N k )e^) 



where (3 S — > 0, so that 

d(N k ) 



Sk = (ftWs) (22) 



dt 

Here, /3 S = l/kT s with T s the temperature of the source, and y^(/3 s ) = al(/3 s )e~^ s k 
where ct k {f3 s ) is the probability per time unit that a phonon is emitted from the mode 
k to the source. 

Finally, by gluing together all the events listed above (points 1. 2. and 3.), we deduce 



d ^ L = s k + M/3) ((N k ) + 1 - (N k )e^) + 

5> W G0) [(N k + l)(iV,.) - (N k )(N 3 + l)e^"^] , k = 1, ...,z (23) 

which correspond exactly to Eqs. (1) with<^(/3) = a fc (/3)e" /3r ^ fc and A jk {P) = A fcj (/3)e /3?lK ^ ) 



C. Qualitative approach to Prohlich condensation 

As mentioned in the Introduction, Frohlich condensation state is a stationary state for 
the system described by Eqs. (??) achieved when the rates of energy supply s k exceed some 
threshold value. In this situation, the energy of the system is primarily located in the mode 
of lowest frequency, i.e., it can be found that (Ni) ~ (N) = /](Nj). 

3 

Although the existence of this state has been emphasized both analytically [2J [TU] and 
numerically [5] from Eqs. (??), Frohlich condensation may be highlighted more directly on 
the basis of the transition probabilities computed in the previous subsection. In that regard, 
it should be stressed that Frohlich condensation is obtained as a consequence of the energy 
redistribution due to the influence of non-linear interactions between the system and the 

10 



heat bath [TO]. Now, we found in the previous subsection that the probability per time unit 
that one phonon is emitted from the mode j while one phonon is simultaneously absorbed 
in the mode k reads as (see Eqs. (??) and (??)) 

W n {n k + \\Uj - l\n k ;nj) = A kj (f3){n k + l)n,-. (24) 

Conversely, the probability per time unit that one phonon is emitted from the mode 
k while one phonon is absorbed in the mode j is simply W n (n k — 1;% + l\n k ;rij) = 
Aj k ((3)n k (rij + 1). Then, using Eq. (??) A jfc , one gets 

W n (n k - 1; nj + l\n k ; n,) = Ajye^-^M^ + 1). (25) 

In particular, assuming that n k ,rij 3> 1, it appears that the probability that one phonon 
is absorbed in the mode of lower frequency through non-linear events is always greater than 
the probability that one phonon is emitted from that mode. Indeed, when n k ,rij ^> 1, one 
gets (n k + l)rij ~ n k (rij + 1). Thus, the ratio between the transition probabilities (??) and 
(??) can be approximated as e fiH(fJk ~ u '\ If Uk < u j: e muJk '^ ] < 1 so that the probability of 
a simultaneous absorption in the mode k and emission from the mode j will be the larger one 
at any time (Eq. (??)). On the contrary, if co k > Uj, e^^-^j) > 1 and the larger probability 
will be that of a simultaneous absorption in the mode mode j and emission from the mode 
k (Eq. (??)). Here, we see how relevant the condition of detailed balance through Eq. (??) 
is in order to get Frohlich condensation. In particular, a wrong choice of the constants A i: ,- 
could lead to a stationary state characterized by the excitation of modes of higher frequency 
in the spectrum of the system. 

Coming back to the number of phonons in each mode, the condition for large values of 
n k , rij can be fulfilled in so far as the rates of energy supply s k are large enough. More 
specifically, it was shown at the end of the last subsection, that the interactions between the 
system and the external source could be depicted similarly to the linear interactions between 
the system and the heat bath with a condition of infinite temperature. For the latter, let us 
recall from subsection II A that the probability that one phonon is emitted from the mode 



k (through linear events) reads as 



W n (n k - l\n k ) = a k (/3)n k , (26) 

11 



while the probability of absorption is given by 

W n (n k + IK) = a k ((3)(n k + \)e~^K (27) 

Here, in the absence of source, the condition for large values for the n k s can never be 
fulfilled because if the n k s become too large, emission will be favored against absorption. 
On the contrary, when only the external source is considered, one has (3 = (3 S — > so that 
according to (??) and (??) (here a k must be changed in s k ) absorption is always favored. 
Thus, in considering the presence of both heat bath and external source, the possibility of 
working with large number of phonons in each normal mode will depend to some extent on 
whether the ratio s k /a k is large or not. In this way, supposing that s k ^> a k for all k, one 
will naturally get n k ^> 1 after a certain time so that non-linear interactions will redistribute 
the energy into the mode of lowest frequency. 

III. THE ORIGINAL FROHLICH MASTER EQUATION AND THE RELATION 
TO THE WU-AUSTIN HAMILTONIAN DESCRIPTION 

In the previous section, we have been able to derive the Frohlich rate equations (??) on 
the basis of purely semi-classical arguments. To summarize, it was assumed that : 

1. The number of phonons in each normal mode of the system (as a function of time) 
could be represented as a (homogeneous) Markov process. This property allowed us to 
describe the evolution of the system on the basis of a classical master equation (??). 

2. Due to the discrete nature of the energy (phonons), each process could be described 
as a birth and death process. 

3. On each normal mode, all phonons could be considered independent. In addition, 
according to the detailed balance condition and the requirement of Bose-Einstein sta- 
tionary distributions of phonons when only one of both the external source and the 
heat bath is present, the transition probabilities that take place in the master equation 
have been worked out. 

In the end, the previous calculations allowed us to derive not only the Frohlich rate equations 
but also the original classical master equation from which the former can be deduced 
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d t p(n, t) = d t p s (n, t) + d t p a (n, t) + d t p A (n, t). (28) 

Here, p s (n, t), p a (n, t) and Pa(ti, t) are the probability density functions of the number of 
phonons in all normal modes of the system (let us emphasize again that n = (ni, n%, ..., n z )) 
whose evolutions are respectively due to : 

• (Linear) events resulting from interactions between the system and the external source; 

• Linear events (i) resulting from interactions between the system and the heat bath; 

• Non-linear events (ii) resulting from interactions between the system and the heat 
bath. 

By substituting Eqs. (??) and (??) into Eq. (??), we get the master equation for p a (n,t) 

d t p a (n,t) = y^ otiip) {[{ni + l)pn{rii + 1, t) - njp n {nj, t)\ + 

(29) 

+e-^[n i p n {n l - l,t) - (m + l)p n (n h t)]} . 
The master equation for p s (n, t) is found in a similar way by letting j3 = j3 s — > and 
changing CKj in Sj 

d t p s {n,t) = ^ Sj(/3 S ) {{rii + l)p n (ni + l,t) + Tiip n (ni - l,t) - (2ra + l)p n {n h t)} . (30) 



Let us recall from subsection II B that p n {rrii } t) is the probability of getting n\ phonons 
in the mode of frequency wi, ri2 phonons in the mode of frequency ix>2,..., n z phonons in the 
mode of frequency u z but rrii phonons (instead of rti phonons) in the mode of frequency u^. 
A similar notation p n {mi] rrij,t) is used when the exception involves two modes. Naturally, 
Pn( n ii t) and p n {rii\ rij,t) correspond to p(n, t). 

Finally, the master equation corresponding to non-linear events is given by substituting 
Eqs. (??) and (??) into Eq. (??) so that 

d t p A (n, t) = ^ Aa(P)ni(nj + l)p n {rii - 1; n 3 -, + 1, t) - 

Aij(p)(rii + l)n j p n (n i ;nj,i) + 

(31) 

Aji(J3)(m + l)njPn(rii + 1; rij - 1, t) - 
Aji(p)rii(nj + l)p n (m;nj,t), 
13 



Eq. (??) together with Eqs. (??), (??) and (??) are the original classical master equations 
from which the Frohlich rate equations (??) may be deduced. In the Introduction, it was 
mentioned that Eqs (??) could be derived from a quantum Hamiltonian that was originally 
introduced by Wu and Austin [U [5] . More specifically, it can be found that the above master 
equations are no more than an intermediate step of that derivation. In that regard, let us 
also emphasized that the possibility of deducing classical master equations (also called Pauli 
master equations) from a microscopic Hamiltonian is well-known in the realm of quantum 
open systems [151 Q5] . In the case of Frohlich theory, that correspondence has not been given 
much attention in the literature. From this point of view as well as to provide a consistency 
check of Eqs. (??) to (??), let us emphasize how the Frohlich master equations can be 
deduced from the quantum approach given by Wu and Austin. 

To that purpose, we recall that the Wu-Austin Hamiltonian [I] is given by 



H = Hn + H,, 



inti 



where H stands for the Hamiltonian of the system of normal modes, the heat bath plus 
the external source in the absence of interactions, and has the following form : 



■-b 



Ho = 2^ ftuia\ a i + ^2 ^kb{h + 22 ftQ'p c l c p- 
i=i fc=i p=i 

Here, a], Oj are the creation and annihilation operators associated with the normal mode 
of frequency Wj of the system, b k ,bk are the creation and annihilation operators associated 
with the normal mode of frequency f^ of the heat bath, and cL c p are the creation and 
annihilation operator associated with the normal mode of frequency Q' of the source. In 
parallel, Hi nt is the interaction Hamiltonian such that 

H^t = ^2 ^hkalh + 22 Xi,j,k a i a )bp + X^ &p a l c p + H - c -> ( 32 ) 

i,k i,j,k *iP 

Since one is interested in the statistical properties of observables related to the system 
only, it is convenient to work with the reduced density operator S(t), which is defined as 
the usual density operator traced over the states related to the heat bath and the source 
[15]. Assuming that the dynamics of each normal mode verifies the Markov property, it may 
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be shown JT5J [T7] that second perturbation theory in the interaction picture applied to the 

interaction Hamiltonian (??) actually yields 



dS 
~dt 



- - ^ { (^ r * + 0i ) \ a i a \ s + Sa i a l ~ SalSaA + 

i 

[0i( r i + 1) + 0"i] ( alai'S' + S'ajaj — 2aiSa\ j > 



9 L^t h2XhJ,J-iXi>,j',j'-i> X 



i,j,i',j' 
j>i,j'>i' 






f ajaJa. ; /a],,S' + Saia\jajia\, — 2ajial,Sa i al j J (r\,-_j + 1) + 
f ajialidiajS + Saj/a^aiOj — 2aia]jSajia i , J r.,_ 



where 



0* = £ (lAvl 2 + |£m| 2 ) , * = ^I6,i| 2 (<<&>. - (6j6,>») • (34) 

Here, (...) s and (...)& are averages performed over states related to the source and the heat 
bath respectively. Now, the probability of observing m phonons in the first mode, ri2 phonons 
in the second, and so on, corresponds to the matrix element (n\S(t)\n) := p(n,t), where 
\n) = \ni)...\n z ). Naturally, |n») is the eigenvector of the operator hi = a^o*. Considering 
only the first sum in Eq. (??) (that stems from the linear terms in the Hamiltonian (??)), 
one can easily deduce the equation : 



d t p(n,t) = - y^(0jrj + o-j) [{uj + J)p n (ni, t) - njp n (nj - 1, t)\ + 

i (35) 

+(<pi(ri + 1) + en) [nip n (ni,t) - (m + l)p n {rii + l,t)] 
In particular when no source is present, one simply gets 

71 | 71 1 71 e P^i 

6iTi + <j; = — |Aj i\ 2 (b]bi) b = — -|Aj A 2 — wz and ddr* + 1) + o; = — -|Aj A 2 — jt , 

as the heat bath is required to be at thermal equilibrium with a temperature /3 = 1/kT. 
Using these two last expressions into Eq. (??), we obtain : 

^l^.^l 2 ^^ -{[{ni + l)p n {rii + l,t) -nip n (m,t)] + 

i h eP**-\ ( 36 ) 

+e- phWl [n i p n (n i -l,t)-(m+ l)p n (m,t)]} 
15 



which is exactly the classical master equation (??) with (Xi(j3) = — | A ^ j | 2 — 37 . Nat- 

ft c- X 

urally, the classical master equation (??) related to the source is found in a similar way by 

considering only the £-terms in Eq. (??) and by using the condition (c{q) s = — 3-r 

with [5 s hjJi <C 1. 

Although no further assumption has been needed to get Eqs. (??) and (??) from Eq. 
(??), the master equation (??) related to non-linear interactions may only be deduced by 
first supposing that all quantum coherences have been destroyed in the system, i.e : 



(n\S(t)\m) = p(n,t)5 ntTn , Vn,m (37) 

where 5 n ^ m = S nitmi 5 n2 ^ m2 ...5 n ^ mz . Surprisingly although Eq. (??) is a quite usual 
assumption when deriving classical master equations from a microscopic Hamiltonian, it 
was never clearly put in evidence in the case of Frohlich theory. 

Finally, one gets (we recall that according to Eq. (??), the condition i — 1 = j — j' is 
required) : 

{n\aia^ajia\,S{t)\n) = 5 i:i i5jj>(ni + l)njp n (ni;nj,t) = (n|S'(£)a i ajaj'aj / |n), 



{n\ajia\,Saia}-\n) = Si^Sjjin^rij + l)p n (ni — l;rij + l,t), 



{n\ajial,aiOjS(t)\n) = Si^djj/n^rij + \)p n {rii\ Uj , t) = (n|S'(t)a : ,-/aJ,aja^|n), 
and 

(n\aia!jSaj>al,\n) = 5 it i>5j tj >(ni + i)rijp n (ni + l;rij - l,t). 
From the second sum in Eq. (??) (terms in x) an d the relations 7\,_, 



and r-j-i + 1 = — ^ r , one obtains : 
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7T 1 

dtp(n,t) ---- J2 ^IXijj-il 2 ^^) _ i {{rii + l)n jPn (ni; n^t) + 

(38) 



riiirij + l)p n {rii - 1; rij + 1, t) + 



efihCi-ui) (( n . + i)n jPn (m + 1; n, - 1, t) 

n i (n j + l)p n {ni, n j} t))}, 
that is exactly the master equation (??) accounting for the non-linear coupling between the 

7T 1 7T e PHui-Uj) 

normal modes with Aj,(/3) = -ttIyj ,• j-i\ 2 — ^-, > andA,j(/3) = t^tIy* ,• ,-_il 2 — ; 5^ ; 

Therefore, we see that the condition for detailed balance [Eq. (??)] is adequately verified. 
As a final complement, let us mention that the Frohlich master equation can be deduced 
not only from the Hamiltonian (??) but also from a large family of Hamiltonians for which 
the operators describing the thermal bath and the source can be very complex objects as 
various combinations of bosonic and/or fermionic, creation and/or annihilation operators 
[IB"] . In the end, only the expressions of the parameters s, if and \ will differ from the ones 
given in this section. We refer the reader interested in a further detailed investigation on 
the quantum Hamiltonian formulation of Frohlich theory to references jlHZ] ■ 



IV. NUMERICAL COMPUTATIONS 

It should be stressed that most theoretical studies on Frohlich condensation are based on 
the numerical integration of the rate equations (??) whereas the issue of the fluctuations of 
the number of phonons in each mode is largely unexplored. From this point of view, a study 
of the distributions of the occupation numbers could be of utmost importance to identify 
Frohlich condensation experimentally at least in qualitative terms. As emphasized above, 
the time evolution of these distributions is provided by Eqs. (??) to (??). These equations, 
being master equations with many variables, are extremely difficult to solve analytically as 
the number of involved processes (here the number of normal modes) increases. Nevertheless, 
several algorithms allow to estimate the solution of such equations numerically. The most 
popular one, proposed by Gillespie, consists in generating stochastic trajectories governed 
by the master equation on the basis of Monte-Carlo procedures [19] . Eventually, these 
trajectories can be used to estimate the expected distribution at any time t as well as the 
related moments. 
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In this context, we investigate numerically Frohlich condensation through two different 
approaches. The first one, as a consistency check, involves the integration of the rate 
equations (??) so that the evolution of the average number of phonons of each normal mode 
is followed for different values of the source parameters Sk- The second approach consists 
in generating stochastic trajectories whose evolution is governed by the Frohlich master 
equation [Eqs. (??) to (??)], by means of the Gillespie algorithm. By collecting the results 
in the form of histograms, it is thus possible to estimate the distribution of phonons for any 
normal mode at any time t. 

Some representative results of Frohlich condensation using the first approach are reported 
in Figure [TJ The system is initially at thermal equilibrium, i.e, the average number of 
phonons is given at t = by the Planck formula : 

W = °)> = ^T7> Vfc - (39) 

Then, each solution tends to a stationary state with a convergence rate that is mainly 
depending on the values of the coupling constants Sk, ctki Ajy. Since we are interested 
in a qualitative description of the condensation, we are not committed to an exhaustive 
exploration of the parameter space (in particular, the coupling constants are given in (time 
arbitrary units) -1 ) [20]. Let us simply specify that we required Ay <C a^ Vfc, j because events 
related to linear interactions are more likely to arise than events related to non-linear ones. 
Also, the coefficients Ay have been fixed according to the relation A^. = Aye^ ^ fe-£J -^ which 



follows a theoretical condition of detailed balance (see section II B ). Now, we see from Figure 
[I] how Frohlich condensation takes place in average terms : a situation wherein Sk < cik, for 
all k, leads to a stationary state relatively close to the initial (thermal equilibrium) condition 
(see Figure flTa)). Alternatively, a situation where Sk 3> ctk holds for the majority of normal 
modes leads to a stationary state characterized by a large number of phonons populating 
the mode of lowest frequency, i.e., (iVi(oo)) ^> (iV&(oo)) for all ft > 1 (see Figure [I] (b)). 
Besides, the values of (iVfc(oo)) with k > 1 still remain close to the thermal equilibrium 
values. 

Representative results of Frohlich condensation using the second approach are depicted 
in Figure 2, 3 and 4. Each figure (histogram) was obtained by generating 10000 stochastic 
trajectories and by storing the values of the number of phonons of a particular mode at a 
time t when the system is considered to have reached the stationary state. In that regard, t 
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FIG. 1. (Color online) Time evolutions of the average number of phonons (Nk(t)) in the three 
modes of lowest frequency of a system of four modes obtained by integrating Eqs. (??)• The 
frequencies of the normal modes are Wj = i.10 Hz |22j . i = 1, ..., 4 whereas the coupling constants 
associated with the interactions between the system and the heat bath are given by a^ = 10 for all 
k and Ajy = 1 for k < j (let us recall that the relation Aj^ = Ajye^ fi ' a,fe— ^'J is required). Initially, 
the system is supposed to be at thermal equilibrium so that {Nk(t = 0)) is given by Planck formula 
(??). The main change between both figures is the rate of energy supply : (a) s& = 1, (b) Sk = 
20, for all k. All coupling constants a, A and s are given in (time arbitrary unit) - , as discussed 
in the text. 



19 



0.01 



0.008 



:k 0.006 

CO 
-Q 

A" 0.004 



0.002 







(a) 



100 200 300 400 500 600 700 

Number of phonons 




50 100 150 200 250 300 350 400 

Number of phonons 



FIG. 2. (Color online) Normalized histograms of the number of phonons in the first two modes of a 
system of four modes at t = 30 arbitrary units (= stationary state) : (a) mode of frequency u>\ (b) 
mode of frequency 002- Both figures have been obtained by generating 10000 stochastic trajectories 
by means of the Gillespie algorithm applied on Eqs. (??) to (??). Initially, the system is supposed 
to be at thermal equilibrium (see text). The frequency of the normal modes as well as the coupling 
constants related to the system-bath interactions are the same as those used Figure 1. The rate of 
energy supply is s^ = 5 (time arbitrary unit) - for all k (for (a) and (b)). 
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FIG. 3. (Color online) Same as Figure 2 with a rate of energy supply Sk = 20 (time arbitrary 
unit)" for all k (for (a) and (b)). 



has been approximated beforehand using the first approach. Again, the system is supposed 
to be initially at thermal equilibrium which means that each initial condition is chosen at 
random following a Bose-Einstein distribution for all the modes 
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FIG. 4. (Color online) Same as Figure 2 with a rate of energy supply s k = 50 (time arbitrary 
unit) - for all k (for (a) and (b)). 



p(n k ,t = 0) = e~^ knk , with Z k = -=—, 

Fy ' Z k 1 - e-P r ^ 

for k = 1, ...,z (rejection sampling). From a general point of view, the averages of the 
number of phonons deduced from the stochastic trajectories are in excellent agreement with 
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those stemming from the integration of the rate equations (??). To this respect, let us recall 
that the rate equations (??) can be deduced from the Frohlich master equation only if the 
correlations between the normal modes can be neglected p(n,t) ~ 1 lp(nj,t) (see section 



II B). So far, with our choice of parameters, this assumption is thus confirmed independently 
of whether Frohlich condensation is reached or not. 

Now, we see from Figure 2, 3 and 4 how Frohlich condensation appears statistically : 
for weak values of Sk (Figure 2) the distributions of the number of phonons in the first two 
modes of lowest frequency (four modes overall, see Figures) remain close to the Bose-Einstein 
distributions. Progressively, by increasing Sk, we see in the first mode that not only is the 
average number of phonons increasing as expected but also does the variance. As a general 
comment concerning the numerical results displayed here, let us remark that there is an 
analogy with the change of the statistical distributions of photon occupation numbers below 
and above the laser transition. In fact, the photon statistics below the laser transition is a 
thermal equilibrium one peaked at zero, whereas, above the transition, the photon statistics 
gives a non-zero most probable value of the number of photons in the lasing mode with a well- 
known Poisson distribution. Loosely speaking something similar occurs in correspondence 
with the condensation transition here explored. Here, we are not interested in deepening 
this analogy beyond the simple remark that the Frohlich condensation transition displays a 
phenomenology typical of non-equilibrium phase transitions in open systems. 

V. CONCLUSIONS 

In the present work, the phenomenon of Frohlich condensation has been addressed in a 
semi-classical framework. With "semi-classical" it is meant that the evolution of the system 
is described by means of classical equations with the addition of energy quantization. In the 
framework of open systems, this leads to the classical master equation replacing the evolution 
equation for the quantum density matrix. The quantum approaches hitherto proposed are 
based on the microscopic Hamiltonian originally proposed by Wu and Austin or on slightly 
modified versions of it. Here, we have proposed a different approach to highlight what are 
the necessary hypotheses that are required to yield the Frohlich condensation. This is not an 
academic exercise because Frohlich phenomenon could be relevant in very different physical 
contexts, and, in particular, to biologically relevant systems undergoing self-organization at 
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different scales (from single macromolecules up to collections of cells). Thus, the approach 
proposed here on the one hand has the merit of a critical reviewing of the physics underlying 
Frohlich condensation and on the other hand has the advantage of providing the time evo- 
lution of the probability density function which allows the computation of higher moments 
of the number of phonons beside their average quantities. Actually, though for a system 
with a small number of normal modes, we have shown a clearcut change of the histograms 
of the number of phonons in correspondence with the condensation transition. Though in a 
somewhat loose way, an analogy with the change of the distribution of photons below and 
above the laser transition has been noted. 

ACKNOWLEDGMENTS 

I warmly thank my colleagues and friends Marco Pettini and Ricardo Lima for enlight- 
ening comments and discussions. 

Appendix A: Detailed Balance and birth and death process 



We consider the homogeneous Markov process Nk(t) of section II A According to the 
birth and death nature of N k (t), the stationary solution p s (n k ) of equation (??) is given 
V rik > by 

= J{n k ) - J{n k - 1), (Al) 

where J{rik) = W{n k \n k + l)p s {n k + 1) — W(rik + l\rik)p s (rik) ■ We now sum (??) so 

o = J2 1 J K) - J fa ~ l )\ = J K) - J (-!)' ( A2 ) 

71j=0 

According conditions (??), J(— 1) = 0, so that the detailed balance condition is obtained 

= J(n k ) = W(n k \n k + l)p s (n k + 1) - W(n k + l\n k )p s (n k ) } Vn t >0 (A3) 
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